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MULTI- ADAPTIVE GALERKIN METHODS FOR ODES II: 
IMPLEMENTATION AND APPLICATIONS* 
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Abstract. Continuing the discussion of the multi-adaptive Galerkin methods mcG{q) and 
mdG((j) presented in [A. Logg, SIAM J. Sci. Comput., 24 (2003), pp. 1879-1902], we present 
adaptive algorithms for global error control, iterative solution methods for the discrete equations, 
features of the implementation Tanganyika, and computational results for a variety of ODEs. Ex- 
amples include the Lorenz system, the solar system, and a number of time-dependent PDEs. 
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1. Introduction. In this paper we apply the multi-adaptive Galerkin methods 
mcG(g) and mdG(g), presented in [lU], to a variety of problems chosen to illustrate 
the potential of multi-adaptivity. Throughout this paper, we solve the ODE initial 
value problem 



(1.1) 



u{t) = f{u{t),t), (o,r], 

u(0) = Wo, 



where u : [0,T] M^, / : x (0,T] ^- is a given bounded function that is 
Lipschitz continuous in u, mq G is a given initial condition, and T > a given 
final time. 

We refer to [TU] for a detailed description of the multi-adaptive methods. Here 
we recall that each component Ui (t) of the approximate solution U (t) is a piecewise 
polynomial of degree qi = qi{t) on a partition of (0, T] into Mi subintervals of lengths 
kij = tij — tij-i, j = 1, . . . , Mi. On the interval /.y = [tij-i, Uj], component Ui(t) is 
thus a polynomial of degree qij . 

Before presenting the examples, we discuss adaptive algorithms for global error 
control and iterative solution methods for the discrete equations. We also give a short 
description of the implementation Tanganyika. 

2. Adaptivity. In this section we describe how to use the a posteriori error 
estimates presented in 10 in an adaptive algorithm. 

2.1. A strategy for adaptive error controL The goal of the algorithm is to 
produce an approximate solution U{t) to (jl.ip within a given tolerance TOL for the 
error e{t) = U{t) — u{t) in a given norm || • ||. The adaptive algorithm is based on the 
a posteriori error estimates, presented in TD , of the form 

N Mi 

(2.1) l|e||<EE^S'^V.. 



*Received by the editors May 23, 2001; accepted for publication (in revised form) May 1, 2003; 
published electronically December 5, 2003. 

http: / /www. siam.org/journals/sisc/25-4/38973.html 

t Department of Computational Mathematics, Chalmers University of Technology, SE-412 96 
Goteborg, Sweden (logg@math.chalmers.sc). 



1119 



1120 



ANDERS LOGG 



or 

N 

(2.2) llell <^5.maxfc^'V,,, 

1=1 ^ 

where {sij} are stability weights, {Si} are stability factors (including interpolation 
constants) , local measure of the residual Ri{U, ■) = Ui — f{U, •) of the approx- 

imate solution U{t), and where we have pij = qij for mcG{q) and pij — qij + 1 for 
mdG(q). 

We use (|2.2p to determine the individual time-steps, which should then be chosen 



as 



(2.3) k 



■I] 



/T0L/7V\ 
V S, r,j ) 



We use (|2.ip to evaluate the resulting error at the end of the computation, noting 
that (P?T|) is sharper than (1^ . 

The adaptive algorithm may then be expressed as follows: Given a tolerance 
TOL > 0, make a preliminary guess for the stability factors and then 

(i) solve the primal problem with time-steps based on ()2.3p . 

(ii) solve the dual problem and compute stability factors and stability weights. 

(iii) compute an error estimate E based on (12.11) . 

(iv) ii E < TOL, then stop, and if not, go back to (i). 

Although this seems simple enough, there are some difficulties involved. For 
one thing, choosing the time-steps based on p.Sp may be difficult, since the residual 
depends implicitly on the time-step. Furthermore, we have to choose the proper data 
for the dual problem to obtain a meaningful error estimate. We now discuss these 
issues. 

2.2. Regulating the time-step. To avoid the implicit dependence on fc^ for 
we may try replacing (j2.3p with 



(2.4) % = 



/ TOL/A^ \ 

\Si n^j^i) 



l/pi: 



Following this strategy, if the time-step on an interval is small (and thus also is the 
residual), the time-step for the next interval will be large, so that (|2.4p introduces 
unwanted oscillations in the size of the time-step. We therefore try to be a bit more 
conservative when choosing the time-step to obtain a smoother time-step sequence. 
For (j2.4p to work, the time-steps on adjacent intervals need to be approximately 
the same, and so we may think of choosing the new time-step as the (geometric) 
mean value of the previous time-step and the time-step given by ()2.4p . This works 
surprisingly well for many problems, meaning that the resulting time-step sequences 
are comparable to what can be obtained with more advanced regulators. 

We have also used standard PID (or just PI) regulators from control theory with 
the goal of satisfying 

(2.5) = T0L/7V, 

or, taking the logarithm with C,, = log(TOL/(A^S'i)), 



(2.6) 



Pij log ki J + lOgTy = Cj, 
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with maximal time-steps {fcjj}, following work by Soderlind il5j and Gustafsson, 
Lundh, and Soderlind [B]. This type of regulator performs a little better than the 
simple approach described above, provided the parameters of the regulator are well 
tuned. 

2.3. Choosing data for the dual. Different choices of data fx and g for the 
dual problem give different error estimates, as described in [lOj . where estimates for 
the quantity 



were derived. The simplest choices are 5 = and {^pT)i — Sin for control of the final 
time error of the nth. component. For control of the ^^-norm of the error at final 
time, we take g — Q and ^pT = c(r)/|je(T)|j with an approximation e of the error e. 
Another possibility is to take ipr = and gt (t) — 6.^ for control of the average error 
in component n. 

If the data for the dual problem are incorrect, the error estimate may also be 
incorrect: with (pT (or g) orthogonal to the error, the error representation gives only 
< TOL. In practice, however, the dual — or at least the stability factors — seems 
to be quite insensitive to the choice of data for many problems so that it is, in fact, 
possible to guess the data for the dual. 

2.4. Adaptive quadrature. In practice, integrals included in the formulation 
of the two methods mcG{q) and mdG((7) have to be evaluated using numerical quadra- 
ture. To control the resulting quadrature error, the quadrature rule can be chosen 
adaptively, based on estimates of the quadrature error presented in [TU] . 

2.5. Adaptive order, q-adaptivity. The formulations of the methods include 
the possibility of individual and varying orders for the different components, as well 
as different and varying time-steps. The method is thus q-adaptive (or p-adaptive) in 
the sense that the order can be changed. At this stage, however, lacking a strategy 
for when to increase the order and when to decrease the time-step, the polynomial 
orders have to be chosen in an ad hoc fashion for every interval. One way to choose 
time-steps and orders could be to solve over some short time-interval with different 
time-steps and orders, and optimize the choice of time-steps and orders with respect 
to the computational time required for achieving a certain accuracy. If we suspect 
that the problem will change character, we will have to repeat this procedure at a 
number of control points. 

3. Solving the discrete equations. In this section we discuss how to solve the 
discrete equations that we obtain when we discretize (|l.ip with the multi-adaptive 
Galerkin methods. We do this in two steps. First, we present a simple explicit 
strategy, and then we extend this strategy to an iterative method. 

3.1. A simple strategy. As discussed in [10], the nonlinear discrete algebraic 
equations for the mcG((7) method (including numerical quadrature) to be solved on 
every local interval lij take the form 



(3.1) e.j™=e.jO + %E^™" /.(f/(r^'(s^')),r,7i(4^1)), m=l,...,<zy. 



n=0 

where {£,ijm}m=i ^^^^ degrees of freedom to be determined for component Ui{t) 
on interval /y, {wmn'lmLi n=o are weights, {sn'^'}^i'4o ^^'^ quadrature points, and Tij 
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maps lij to (0, 1]: Tij{t) — {t — ti,j-i)/ {tij ^ Uj-i)- The discrete equations for the 
nidG((7) method are similar in structure and so we focus on the mcG^q) method. 

The equations are conveniently written in fixed point form, so we may apply 
fixed point iteration directly to p.l[) : i.e., we make an initial guess for the values of 
{CymlmLii e.g., £,ijm = £,ijo for TO = 1, ... , Qij , and then compute new values for these 
coefficients from p.ip . repeating the procedure until convergence. 

Note that component Ui (t) is coupled to all other components through the right- 
hand side fi — fi{U,-). This means that we have to know the solution for all other 
components in order to compute Ui{t). Conversely, we have to know Ui{t) to compute 
the solutions for all other components, and since all other components step with 
different time-steps, it seems at first very difficult to solve the discrete equations 

m- 

As an initial simple strategy we may try to solve the system of nonlinear equations 
p.ip by direct fixed point iteration. All unknown values, for the component itself and 
all other needed components, are interpolated or extrapolated from their latest known 
values. Thus, if for component i we need to know the value of component I at some 
time ti, and we only know values for component I up to time t; < tt, the strategy is 
to extrapolate Ui (t) from the interval containing ti to time ti , according to the order 
of Ui{t) on that interval. 

In what order should the components now make their steps? Clearly, to update a 
certain component on a specific interval, we would like to use the best possible values 
of the other components. This naturally leads to the following strategy: 

(3.2) The last component steps first. 

This means that we should always make a step with the component that is closest 
to time t = 0. Eventually (or after one step), this component catches up with one 
of the other components, which then in turn becomes the last component, and the 
procedure continues according to the strategy p.2p . as described in Figure [3Jl 

This gives an explicit time-stepping method in which each component is updated 
individually once, following (|3.2p . and in which we never go back to correct mistakes. 
This corresponds to fixed point iterating once on the discrete equations p.ip . which 
implicitly define the solution. We now describe how to extend this explicit time- 
stepping strategy to an iterative process, in which the discrete equations are solved 
to within a prescribed tolerance. 

3.2. An iterative method. To extend the explicit strategy described in the 
previous section to an iterative method, we need to be able to go back and redo 
iterations if necessary. We do this by arranging the elements — we think of an element 
as a component Ui(t) on a local interval hj — in a time-slab. This contains a number of 
elements, a minimum of N elements, and moves forward in time. On every time-slab, 
we have to solve a large system of equations, namely, the system composed of the 
element equations p.ip for every element within the time-slab. We solve this system 
of equations iteratively, by direct fixed point iteration, or by some other method as 
described below, starting from the last element in the time-slab, i.e., the one closest to 
t — 0, and continuing forward to the first element in the time-slab, i.e., the one closest 
to t = T. These iterations are then repeated from beginning to end until convergence, 
which is reached when the computational residuals, as defined in |10j . on all elements 
are small enough. 

3.3. The time-slab. The time-slab can be constructed in many ways. One is by 
dyadic partitioning, in which we compute new time-steps for all components, based 
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T 

Fig. 3.1. The last component steps first and all needed values are extrapolated or interpolated. 



on residuals and stability properties, choose the largest time-step K as the length of 
the new time-slab, and then, for every component, choose the time-step as a fraction 
K/2"'. The advantage of such a partition are that the time-slab has straight edges; i.e., 
for every time-slab there is a common point in time t' (the end-point of the time-slab) 
which is an end-point for the last element of every component in the time-slab, and 
that the components have many common nodes. The disadvantage is that the choice 
of time-steps is constrained. 

Another possibility is a rational partition of the time-slab. We choose the largest 
individual time-step K as the length of the time-slab, and time-steps for the remaining 
components are chosen as fractions of this large time-step, K/2, K/3, K/A, and so on. 
In this way we increase the set of possible time-step selections, as compared to dyadic 
partitioning, but the number of common nodes shared between different components 
is decreased. 

A third option is to not impose any constraint at all on the selection of time- 
steps — except that we match the final time end-point. The time-steps may vary 
also within the time-slab for the individual components. The price we have to pay 
is that we have in general no common nodes, and the edges of the time-slab are 
no longer straight. We illustrate the three choices of partitioning schemes in Figure 
13.21 Although dyadic or rational partitioning is clearly advantageous in terms of easier 
bookkeeping and common nodes, we focus below on unconstrained time-step selection. 
In this way we stay close to the original, general formulation of the multi-adaptive 
methods. We refer to this as the time- crawling approach. 

3.4. The time-crawling approach. The construction of the time-slab brings 
with it a number of technical and algorithmic problems. We will not discuss here the 
implementational and data structural aspects of the algorithm — there will be much 
to keep track of and this has to be done in an efficient way — but we will give a brief 
account of how the time-slab is formed and updated. 
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T 

Fig. 3.2. Different choices of time-slabs. Top: a dyadic partition of the time-slab; middle: a 
rational partition; bottom: a partition used in the time-crawling approach, where the only restriction 
on the time-steps is that we match the final time end-point T. 



Assume that in some way we have formed a time-slab, such as the one in Figure 
13.31 We make iterations on the time-slab, starting with the last element and continuing 
to the right. After iterating through the time-slab a few times, the computational 
(discrete) residuals, corresponding to the solution of the discrete equations p.ip . on 
all elements have decreased below a specified tolerance for the computational error, 
indicating convergence. 

For the elements at the front of the slab (those closest to time t = T), the 
values have been computed using extrapolated values of many of the other elements. 
The strategy now is to leave behind only those elements that are fully covered by all 
other elements. These are then cut off from the time-slab, which then decreases in 
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Fig. 3.3. The time-slab used in the time-crawling approach to multi- adaptive time-stepping 
(dark grey). Light grey indicates elements that have already been computed. 



size. Before starting the iterations again, we have to form a new time-slab. This wiU 
contain the elements ol the old time-slab that were not removed, and a number of new 
elements. We form the new time-slab by requiring that all elements of the previous 
time-slab be totally covered within the new time-slab. In this way we know that every 
new time-slab will produce at least N new elements. The time-slab is thus crawling 
forward in time rather than marching. 

An implementation of the method then contains the three consecutive steps de- 
scribed above: iterating on an existing time-slab, decreasing the size of the time-slab 
(cutting off elements at the end of the time-slab, i.e., those closest to time t — 0), and 
incorporating new elements at the front of the time-slab. 

Remark 3.1. Even if an element within the time-slab is totally covered by all 
other elements, the values on this element still may not be completely determined, if 
they are based on the values of some other element that is not totally covered, or if 
this element is based on yet another element that is not totally covered, and so on. To 
avoid this, one can impose the requirement that the time-slabs should have straight 
edges. 

3.5. Diagonal Nevifton. For stiff problems the time-step condition required for 
convergence of direct fixed point iteration is too restrictive, and we need to use a more 
implicit solution strategy. 

Applying a full Newton's method, we increase the range of allowed time-steps and 
also the convergence rate, but this is costly in terms of memory and computational 
time, which is especially important in our setting, since the size of the slab may often 
be much larger than the number of components, N (see Figure [373)) . We thus look for 
a simplified Newton's method which does not increase the cost of solving the problem, 
as compared to direct fixed point iteration, but still has some advantages of the full 
Newton's method. 

Consider for simplicity the case of the multi-adaptive backward Euler method, 
i.e., the mdG(O) method with end-point quadrature. On every element we then want 
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to solve 

(3.3) Uij = Uij^i + kijfi{U{tij),Uj). 
In order to apply Newton's method we write p.3p as 

(3.4) F{V) = 

with Fi{V) — Uij — Ui.j-i — kij fi{U{tij),tij) and Vi ~ Uij. Newton's method is then 

(3.5) = V"- {F'{V''))-^F{V"). 

We now simply replace the Jacobian with its diagonal so that for component i we 
have 

(3.6) = t/" - ~ 

^ ' ij V 1 _ i. Mi. 

^ dUi 

with the right-hand side evaluated at We now note that we can rewrite this as 

(3.7) - [/- - e{ur, - CA,,-i - h,h) = (1 - 0)ur, + ^(c/.,,-i + h,h) 

with 
(3.8) 



1 ^ij di'i 

so that we may view the simplified Newton's method as a damped version, with 
damping 0, of the original fixed point iteration. 

The individual damping parameters are cheap to compute. We do not need to 
store the Jacobian and we do not need linear algebra. We still obtain some of the 
good properties of the full Newton's method. 

For the general mcG(g) or mdG(g) method, the same analysis applies. In this 
case, however, when we have more degrees of freedom to solve for on every local 
element, 1 — fcy |^ will be a small local matrix of size q x q for the nicG(q) method 
and size {q + 1) x {q + 1) for the mdG(q) method. 

3.6. Explicit or implicit. Both mcG(g) and mdG((7) are implicit methods since 
they are implicitly defined by the set of equations (13. ip on each time-slab. However, 
depending on the solution strategy for these equations, the resulting fully discrete 
scheme may be of more or less explicit character. Using a diagonal Newton's method 
as in the current implementation of Tanganyika, we obtain a method of basically 
explicit nature. This gives an efhcient code for many applications, but we may expect 
to meet difhculties for stiff problems. 

3.7. The stiffness problem. In a stiff problem the solution varies quickly inside 
transients and slowly outside transients. For accuracy the time-steps will be adap- 
tively kept small inside transients and then will be within the stability limits of an 
explicit method, while outside transients larger time-steps will be used. Outside the 
transients the diagonal Newton's method handles stiff problems of sufhciently diago- 
nal nature. Otherwise the strategy is to decrease the time-steps whenever needed for 
stability reasons. Typically this results in an oscillating sequence of time-steps where 
a small number of large time-steps are followed by a small number of stabilizing small 
time-steps. 
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Our solver Tanganyika thus performs like a modern unstable jet fighter, which 
needs small stabilizing wing flaps to follow a smooth trajectory. The pertinent ques- 
tion is then the number of small stabilizing time-steps per large time-step. We ana- 
lyze this question in 0] and show that for certain classes of stiff problems it is indeed 
possible to successfully use a stabilized explicit method of the form implemented in 
Tanganyika. 

3.8. Preparations. There are many "magic numbers" that need to be computed 
in order to implement the multi-adaptive methods, such as quadrature points and 
weights, the polynomial weight functions evaluated at these quadrature points, etc. 
In Tanganyika, these numbers are computed at the startup of the program and stored 
for efficient access. Although somewhat messy to compute, these are all computable 
by standard techniques in numerical analysis; see, e.g., |13) . 

3.9. Solving the dual problem. In addition to solving the primal problem 
(|1.1[) . we also have to solve the continuous dual problem to obtain error control. This 
is an ODE in itself that we can solve using the same solver as for the primal problem. 

In order to solve this ODE, we need to know the Jacobian of / evaluated at a 
mean value of the true solution u{t) and the approximate solution U{t). If U{t) is 
sufhciently close to u, which we will assume, we approximate the (unknown) mean 
value by U{t). When solving the dual, the primal solution must be accessible, and the 
Jacobian must be computed numerically by difference quotients if it is not explicitly 
known. This makes the computation of the dual solution expensive. Error control 
can, however, be obtained at a reasonable cost: for one thing, the dual problem does 
not have to be solved with as high a precision as the primal problem; a relative error 
of, say, 10% may be disastrous for the primal, whereas for the dual this only means 
that the error estimate will be off by 10%, which is acceptable. Second, the dual 
problem is linear, which may be taken into account when implementing a solver for 
the dual. If we can afford the linear algebra, as we can for reasonably small systems, 
we can solve the discrete equations directly without any iterations. 

4. Tanganyika. We now give a short description of the implementation of the 
multi-adaptive methods, Tanganyika, which has been used to obtain the numerical 
results presented below. 

4.1. Purpose. The purpose of Tanganyika [12] is to be a working implementa- 
tion of the multi-adaptive methods. The code is open-source (GNU GPL p!]), which 
means that anyone can freely review the code, which is available at http: / /www.phij 
chalmers.se/tanganyika/. Comments are welcome. 

4.2. Structure and implementation. The solver is implemented as a C/C-I--I- 
library. The CH — h language makes abstraction easy, which allows the implementation 
to follow closely the formulation of the two methods. Different objects and algorithms 
are thus implemented as CH — h classes, including Solution, Element, cGqElement, 
dGqElement, TimeSlab, ErrorControl, Galerkin, Component, and so on. 

5. Applications. In this section, we present numerical results for a variety of 
applications. We discuss some of the problems in detail and give a brief overview of 
the rest. A more extensive account can be found in [llj . 

5.1. A simple test problem. To demonstrate the potential of the multi-adaptive 
methods, we consider a dynamical system in which a small part of the system oscil- 
lates rapidly. The problem is to accurately compute the positions (and velocities) of 
the N point-masses attached with springs of equal stiffness, as in Figure [5TT1 
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Fig. 5.1. A mechanical system consisting of N = 5 masses attached with springs. 



If we choose one of the masses to be much smaller than the others, mi = 10" 
and 7712 — 7713 — ■ ■ ■ — niN ~ then we expect the dynamics of the system to be 
dominated by the smallest mass, in the sense that the resolution needed to compute 
the solution will be completely determined by the fast oscillations of the smallest 
mass. 

To compare the multi-adaptive method with a standard method, we first compute 
with constant time-steps fc = fco using the standard cG(l) method and measure the 
error, the cpu time needed to obtain the solution, the total number of steps, i.e., 
M = J2iLi and the number of local function evaluations. We then compute the 
solution with individual time-steps, using the mcG(l) method, choosing the time-steps 
ki ~ fco for the position and velocity components of the smallest mass, and choosing 
ki = lOOfco for the other components (knowing that the frequency of the oscillations 
scales like 1 / y/rn). For demonstration purposes, we thus choose the time-steps a priori 
to fit the dynamics of the system. 

We repeat the experiment for increasing values of N (see Figure 15. 2p and find 
that the error is about the same and constant for both methods. As N increases, the 
total number of time-steps, the number of local function evaluations (including also 
residual evaluations), and the cpu time increase linearly for the standard method, as 
we expect. For the multi-adaptive method, on the other hand, the total number of 
time-steps and local function evaluations remains practically constant as we increase 
N. The cpu time increases somewhat, since the increasing size of the time-slabs 
introduces some overhead, although not nearly as much as for the standard method. 
For this particular problem the gain of the multi-adaptive method is thus a factor TV, 
where N is the size of the system, so that by considering a large-enough system, the 
gain is arbitrarily large. 



5.2. The Lorenz system. We consider now the famous Lorenz system. 



(5.1) 



y — rx ~ y ~ xz, 
z = xy — bz, 



with the usual data (a;(0), 2/(0), z(0)) = (1,0,0), a = 10, b = 8/3, and r = 28; see 
[5]. The solution u{t) — {x{t),y{t), z{t)) is very sensitive to perturbations and is often 
described as "chaotic." With our perspective this is reflected by stability factors with 
rapid growth in time. 
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Fig. 5.2. Error, cpu time, total number of steps, and number of function evaluations as function 
of the number of masses for the multi- adaptive cG(l) method (dashed lines) and the standard cG(l) 
method (solid lines). 




The computational challenge is to solve the Lorenz system accurately on a time- 
interval [0, T] with T as large as possible. In Figure [531 is shown a computed solution 
which is accurate on the interval [0,40]. We investigate the computability of the 
Lorenz system by solving the dual problem and computing stability factors to find 
the maximum value of T. The focus in this section is not on multi-adaptivity — we 
will use the same time-steps for all components, and so racG{q) becomes cG{q) — but 
on higher-order methods and the precise information that can be obtained about the 
computability of a system from solving the dual problem and computing stability 
factors. 

As an illustration, we present in Figure [5]4] solutions obtained with different meth- 
ods and constant time-step fc = 0.1 for all components. For the lower-order methods, 
cG(5) to cG(ll), it is clear that the error decreases as we increase the order. Starting 
with the cG(12) method, however, the error does not decrease as we increase the 
order. To explain this, we note that in every time-step a small round-off error of size 
~ 10^^^ is introduced if the computation is carried out in double precision arithmetic. 
These errors accumulate at a rate determined by the growth of the stability factor 
for the computational error (see [I^). As we shall see below, this stability factor 
grows exponentially for the Lorenz system and reaches a value of 10^^ at final time 
T = 50, and so at this point the accumulation of round-off errors results in a large 
computational error. 
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Fig. 5.4. Solutions for the x-component of the Lorenz system with methods of different order, 
using a constant time-step k = 0.1. Dotted lines indicate the point beyond which the solution is no 
longer accurate. 
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Fig. 5.5. The stability factor for computational and quadrature errors, as function of time for 
the Lorenz system. 



5.2.1. Stability factors. We now investigate the computation of stability fac- 
tors for tlie Lorenz system. For simplicity we focus on global stability factors, such 
as 

(5.2) S'M(T) = max / \\'p^''\t)\\ dt, 

lkll=i Jo 

where ip{t) is the solution of the dual problem obtained with ip{T) = v (and g — 0). 
Letting ^{t) be the fundamental solution of the dual problem, we have 

max / ||$(«)(i)w|l < / max ||$(«)(t)^;|l = / |l$(«)(t)|l di = 5t«l(T). 
\M=T-Jo Jo Iklhi Jo 

(5.3) 

This gives a bound S^'^^T) for S^'^^{T), which for the Lorenz system turns out to be 
quite sharp and which is simpler to compute since we do not have to compute the 
maximum in (15.21) . 

In Figure [53] we plot the growth of the stability factor for q = 0, corresponding to 
computational and quadrature errors as described in [10]. The stability factor grows 
exponentially with time, but not as fast as indicated by an a priori error estimate. 
An a priori error estimate indicates that the stability factors grow as 



(5.4) S'[«l (T) - y^'e 



AT 
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where ^ is a bound for the Jacobian of the right-hand side for the Lorenz system. A 
simple estimate is ^ = 50, which already at T = 1 gives S'["l(r) ~ 10^^. In view of 
this, we would not be able to compute even to T = 1, and certainly not to T = 50, 
where we have S^^\T) « 10^*^™. The point is that although the stability factors grow 
very rapidly on some occasions, such as near the first flip at T = 18, the growth is not 
monotonic. The stability factors thus effectively grow at a moderate exponential rate. 

5.2.2. Conclusions. To predict the computability of the Lorenz system, we 
estimate the growth rate of the stability factors. A simple approximation of this 
growth rate, obtained from numerically computed solutions of the dual problem, is 

(5.5) ^['l(r)«4.10(«-3)+o.37T 
or just 

(5.6) (T) « 10«+'^/^ 

From the a posteriori error estimates presented in [10| , we find that the computational 
error can be estimated as 

(5.7) ^;c«^^°l(T)max||7^^||, 

[0,T] 

where the computational residual Tip is defined as 

(5.8) 7^f (i) = \U{U,) - U{U^,-i) - [ f,{U, ■) dt] . 



k 



With 16 digits of precision a simple estimate for the computational residual is ^r-lO 
which gives the approximation 



^16 



(5.9) 10^/3 10-16 = 10^/3 . 

mm kij mm kij 

With time-steps kij =0.1 as above we then have 

(5.10) Ec w 10^/^"^^ 

and so already at time T = 45 we have Eq ~ 1 and the solution is no longer accurate. 
We thus conclude by examination of the stability factors that it is difficult to reach 
beyond time T = 50 in double precision arithmetic. (With quadruple precision we 
would be able to reach time T — 100.) 

5.3. The solar system. We now consider the solar system, including the Sun, 
the Moon, and the nine planets, which is a particularly important n-body problem of 
the form 

» ^ Grriimj 

(5.11) miX^ = y 7^{xj-Xi), 

^ \Xj - 3 

where Xi{t) — {x} {t) , xf {t) , Xi {t)) denotes the position of body i at time t, rrn is the 
mass of body i, and G is the gravitational constant. 

As initial conditions we take the values at 00.00 Greenwich mean time on January 
1, 2000, obtained from the United States Naval Observatory [2;, with initial velocities 




obtained by fitting a high-degree polynomial to the values of December 1999. This 
initial data should be correct to five or more digits, which is similar to the available 
precision for the masses of the planets. We normalize length and time to have the space 
coordinates per astronomical unit, AU, which is (approximately) the mean distance 
between the Sun and Earth, the time coordinates per year, and the masses per solar 
mass. With this normalization, the gravitational constant is Att'^. 

5.3.1. Predictability. Investigating the predictability of the solar system, the 
question is how far we can accurately compute the solution, given the precision in 
initial data. In order to predict the accumulation rate of errors, we solve the dual 
problem and compute stability factors. Assuming the initial data is correct to five or 
more digits, we find that the solar system is computable on the order of 500 years. 
Including also the Moon, we cannot compute more than a few years. The dual solution 
grows linearly backward in time (see Figure [5^, and so errors in initial data grow 
linearly with time. This means that for every extra digit of increased precision, we 
can reach 10 times further. 

5.3.2. Computability. To touch briefiy on the fundamental question of the 
computability of the solar system, concerning how far the system is computable with 
correct initial data and correct model, we compute the trajectories for Earth, the 
Moon, and the Sun over a period of 50 years, comparing different methods. Since 
errors in initial data grow linearly, we expect numerical errors as well as stability 
factors to grow quadratically. 
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Fig. 5.7. The growth of the error over 50 years for the Earth-Moon-Sun system as described 
in the text. 



In Figure 15.71 we plot the errors for the 18 components of the solution, com- 
puted for k = 0.001 with cG(l), cG(2), dG(l), and dG(2). This figure contains 
much information. To begin with, we see that the error seems to grow linearly for 
the cG methods. This is in accordance with earlier observations [H [7] for periodic 
Hamiltonian systems, recalling that the (m)cG((7) methods conserve energy [1^. The 
stability factors, however, grow quadratically and thus overestimate the error growth 
for this particular problem. In an attempt to give an intuitive explanation of the 
linear growth, we may think of the error introduced at every time-step by an energy- 
conserving method as a pure phase error, and so at every time-step the Moon is 
pushed slightly forward along its trajectory (with the velocity adjusted accordingly). 
Since a pure phase error does not accumulate but stays constant (for a circular orbit), 
the many small phase errors give a total error that grows linearly with time. 

Examining the solutions obtained with the dG(l) and dG(2) methods, we see 
that the error grows quadratically, as we expect. For the dG(l) solution, the error 
reaches a maximum level of ^ 0.5 for the velocity components of the Moon. The error 
in position for the Moon is much smaller. This means that the Moon is still in orbit 
around Earth, the position of which is still very accurate, but the position relative 
to Earth, and thus also the velocity, is incorrect. The error thus grows quadratically 
until it reaches a limit. This effect is also visible for the error of the cG(l) solution; 
the linear growth flattens out as the error reaches the limit. Notice also that even if 
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Fig. 5.8. The growth of the error over 5 years for the Earth-MoonSun system computed with 
the mcG(2) method, together with the multi-adaptive time-steps. 



the higher-order dG(2) method performs better than the cG(l) method on a short 
time-interval, it will be outrun on a long enough interval by the cG(l) method, which 
has linear accumulation of errors (for this particular problem). 

5.3.3. Multi-adaptive time-steps. Solving with the multi-adaptive method 
mcG(2) (see Figure [STS)) . the error grows quadratically. We saw in [10| that in order 
for the mcG(g) method to conserve energy, we require that corresponding position 
and velocity components use the same time-steps. Computing with different time- 
steps for all components, as here, we thus cannot expect to have linear error growth. 
Keeping kfvi < tol with tol = 10"^" as here, the error grows as 10~'^T^ and we are 
able to reach T ^ 100. Decreasing tol to, say, 10~^*, we could instead reach T ^ 10^. 

We investigated the passing of a large comet close to Earth and the Moon and 
found that the stability factors increase dramatically at the time t' when the comet 
comes very close to the Moon. The conclusion is that if wc want to compute accurately 
to a point beyond t' , we have to be much more careful than if we want to compute 
to a point just before t' . This is not evident from the size of residuals or local errors. 
This is an example of a Hamiltonian system for which the error growth is neither 
linear nor quadratic. 



Flc;. 5.9. A space-time plot of the solution, (above) and tune-steps (below) for tlie propagating 
front problem, with time going to the right. The two parts of the plots represent the components for 
the two species Ai (lower parts) and A2 (upper parts). 



5.4. A propagating front problem. The system of PDEs 



(5.12) 



{ 



■iii — eu" = — M1U2, 

U2 — = U1U2 



on (0, 1) X (0,T] with e — 0.00001, T — 100, and homogeneous Neumann boundary 
conditions at a; = and x — 1 models isothermal autocatalytic reactions (see |14) ^ 
Ai + 2A2 A2 + 2A2. We choose the initial conditions as 



with xo = 0.2 and U2{x, 0) = 1 — mi(x, 0). An initial reaction where substance Ai is 
consumed and substance A2 is formed will then occur at x = xq, resulting in a decrease 
in the concentration ui and an increase in the concentration U2. The reaction then 
propagates to the right until all of substance Ai is consumed and we have ui = and 
7X2 = 1 in the entire domain. 

Solving with the mcG(2) method, we find that the time-steps are small only close 
to the reaction front; see Figures [531 and [5.101 The reaction front propagates to the 
right as does the domain of small time-steps. 

It is clear that if the reaction front is localized in space and the domain is large, 
there is a lot to gain by using small time-steps only in this area. To verify this, we 




x < Xq, 
X > Xq, 
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Fig. 5.10. The concentrations of the two species U\ and U2 at time t = 50 as function of space 
(above) and the corresponding time-steps (below). 



compute the solution to within an accuracy of 10^^ for the final time error with con- 
stant time-steps ki{t) = kg for all components and compare with the multi-adaptive 
solution. Computing on a space grid consisting of 16 nodes on [0, 1] (resulting in a 
system of ODEs with 32 components), the solution is computed in 2.1 seconds on a 
regular workstation. Computing on the same grid with the multi-adaptive method 
(to within the same accuracy), we find that the solution is computed in 3.0 seconds. 
More work is thus required to compute the multi-adaptive solution, and the reason 
is the overhead resulting from additional bookkeeping and interpolation in the multi- 
adaptive computation. However, increasing the size of the domain to 32 nodes on 
[0, 2] and keeping the same parameters otherwise, we find the solution is now more 
localized in the domain and we expect the multi-adaptive method to perform better 
compared to a standard method. Indeed, the computation using equal time-steps now 
takes 5.7 seconds, whereas the multi-adaptive solution is computed in 3.4 seconds. In 
the same way as previously shown in section 15.11 adding extra degrees of freedom 
does not substantially increase the cost of solving the problem, since the main work 
is done time-stepping the components, which use small time-steps. 

5.5. Burger's equation with moving nodes. As a final example, we present a 
computation in which we combine multi-adaptivity with the possibility of moving the 
nodes in a space discretization of a time-dependent PDE. Solving Burger's equation. 



(5.13) 



it + /imm' — eu" =^ 0, 
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Fig. 5.11. The solution to Burger's equation as function of space at t = 0, t = T/2, and t = T. 



on (0, 1) X (0,T] with initial condition 



(5.14) 



Uo{x) 



sin(7ra;/xo), < x < Xq, 
elsewhere, 



and with fi — 0.1, e = 0.001, and xq — 0.3, we find the solution is a shock forming near 
X = xo- Allowing individual time-steps within the domain, and moving the nodes of 
the space discretization in the direction of the convection, (1, fiu), we make the ansatz 



N 



(5.15) 



where the {Ci}t=i the individual components computed with the multi-adaptive 
method, and the {'Pi{-,t)}^i are piecewise linear basis functions in space for any 
fixed t. 

Solving with the mdG(O) method, the nodes move into the shock, in the direction 
of the convection, so what we are really solving is a heat equation with multi-adaptive 
time-steps along the streamlines; see Figures [5.111 and 15. 121 

6. Future work. Together with the companion paper [10] (see also Oil]), this 
paper serves as a starting point for further investigation of the multi-adaptive Galerkin 
methods and their properties. 
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Fig. 5.12. Node paths for the multi-adaptive moving-nodes solution of Burger's equation. 

Future work will include a more thorough investigation of the application of 
the multi-adaptive methods to stiff ODEs, as well as the construction of efhcient 
multi-adaptive solvers for time-dependent PDEs, for which memory usage becomes an 
important issue. 
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